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ABSTRACT 


A material model which incorporates several key capabilities which have been 
identified by the aerospace community as lacking in state-of-the art composite impact 
models is under development. In particular, a next generation composite impact 
material model, jointly developed by the FAA and NASA, is being implemented into 
the commercial transient dynamic finite element code LS-DYNA. The material model, 
which incorporates plasticity, damage, and failure, utilizes experimentally based 
tabulated input to define the evolution of plasticity and damage and the initiation of 
failure as opposed to specifying discrete input parameters (such as modulus and 
strength). The plasticity portion of the orthotropic, three-dimensional, macroscopic 
composite constitutive model is based on an extension of the Tsai-Wu composite 
failure model into a generalized yield function with a non-associative flow rule. For 
the damage model, a strain equivalent formulation is utilized to allow for the 
uncoupling of the deformation and damage analyses. In the damage model, a semi- 
coupled approach is employed where the overall damage in a particular coordinate 
direction is assumed to be a multiplicative combination of the damage in that direction 
resulting from the applied loads in the various coordinate directions. Due to the fact 
that the plasticity and damage models are uncoupled, test procedures and methods to 
both characterize the damage model and to covert the material stress-strain curves 
from the true (damaged) stress space to the effective (undamaged) stress space have 
been developed. A methodology has been developed to input the experimentally 
determined composite failure surface in a tabulated manner. An analytical approach is 
then utilized to track how close the current stress state is to the failure surface. 
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INTRODUCTION 


As composite materials are gaining increased use in aircraft components where 
impact resistance 1s critical (such as the turbine engine fan case), the need for accurate 
material models to simulate the deformation, damage and failure response of polymer 
matrix composites under impact conditions is gaining in importance. While there are 
several material models currently available within commercial transient dynamic 
codes such as LS-DYNA [1] to analyze the impact response of composites, areas have 
been identified where the predictive capabilities of these models can be improved. 
Most importantly, the existing models often require extensive correlation based on 
structural level impact tests, which significantly limits the ability to use these models 
as predictive tools. Furthermore, most of the existing models apply either a plasticity 
based approach (such as that used by Sun and Chen [2]) or a continuum damage 
mechanics approach (such as that used by Matzenmiller et al [3]) to simulate the 
nonlinearity that takes place in the composite response. As documented in detail by 
Goldberg, et al [4, 5], either of these approaches can capture certain aspects of the 
actual composite behavior. However, optimally, combining a plasticity based 
deformation model (to capture the rate dependence and significant nonlinearity, 
particularly in shear, observed in the composite response) with a damage model (to 
account for the nonlinear unloading and strain softening observed after the peak stress 
is reached) can provide advantages over using one approach or the other. In addition, 
the input to current material models generally consists of point-wise properties (such 
as the modulus, failure stress or failure strain in a particular coordinate direction) that 
leads to curve fit approximations to the material stress-strain curves. This type of 
approach either leads to models with only a few parameters, which provide a crude 
approximation at best to that actual material response, or to models with many 
parameters which require a large number of complex tests to characterize. An 
improved approach is to use tabulated data, obtained from well-defined set of 
straightforward experiments. Using tabulated data allows the actual material response 
data to be entered in discretized form, which permits a more accurate representation of 
the actual material response. 

To begin to address these needs, a new composite material model is being 
developed and implemented for use within the commercial transient dynamic finite 
element code LS-DYNA. The material model is meant to be a fully generalized 
model suitable for use with any composite architecture (laminated or textile). The 
deformation model is based on extending the commonly used Tsai-Wu composite 
failure model [6] to a strain hardening plasticity model with a non-associative flow 
rule. For the damage model, a strain equivalent formulation is used in which the 
deformation and damage calculations can be uncoupled. A significant feature in the 
developed damage model is that a semi-coupled approach has been utilized in which a 
load in a particular coordinate direction results in damage (and thus stiffness 
reduction) in multiple coordinate directions. While different from the approach used 
in many existing damage mechanics models [3] in which a load in a particular 
coordinate direction only leads to a stiffness reduction in the load direction, this 
approach has the potential to more accurately reflect the damage behavior that actually 
takes place, particularly for composites with more complex fiber architectures. To 
model the failure response of the composite, the stress (or strain) based three 
dimensional failure envelope for the composite is entered in a tabulated fashion. A set 


of stress (or strain) invariants and stress (or strain) ratios are then used to determine if 
the current stress (or strain) state is inside or outside of the failure envelope. 

In the following sections of this paper, a brief summary of the plasticity based 
deformation model is presented. The key aspects of the damage model are then 
described, along with details of the procedures that are required to characterize the 
damage model. The fundamental aspects of the failure model along with a set of 
proposed procedures to characterize the failure model are then discussed. 


DEFORMATION MODEL 


A complete description of the deformation model is given in Goldberg et al [4, 5]. 
A summary of the key features of the model is presented here. In the deformation 
model, a general quadratic three-dimensional orthotropic yield function based on the 
Tsai-Wu failure model is specified as follows, where 1, 2, and 3 refer to the principal 
material directions 
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In the yield function, oj represents the stresses and Fi, and Fx are coefficients that 
vary based on the current values of the yield stresses in the various coordinate 
directions. By allowing the coefficients to vary, the yield surface evolution and 
hardening in each of the material directions can be precisely defined. The values of 
the normal and shear coefficients can be determined by simplifying the yield function 
for the case of unidirectional tensile and compressive loading in each of the coordinate 
directions along with shear tests in each of the shear directions. In the above equation, 
the stresses are the current value of the yield stresses in the normal and shear 
directions. To determine the values of the off-axis coefficients (which are required to 
capture the stress interaction effects), the results from 45° off-axis tests in the various 
coordinate directions can be used. The values of the off-diagonal terms in the yield 
function can also be modified as required in order to ensure that the yield surface is 
convex [5]. 

A non-associative flow rule is used to compute the evolution of the components of 
plastic strain. The plastic potential for the flow rule is shown below 
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where oij are the current values of the stresses and Hi; are independent coefficients, 
which are assumed to remain constant. The values of the coefficients are computed 
based on average plastic Poisson’s ratios [4, 5]. The plastic potential function in 
Equation (2) is used in a flow law to compute the components of the plastic strain rate, 
where the usual normality hypothesis from classical plasticity [7] is assumed to apply 


and the variable 4 is a scalar plastic multiplier. 
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Given the flow law, the principal of the equivalence of plastic work [7] can be 
used to determine that the plastic potential function h can be defined as the effective 
stress and the plastic multiplier can be defined as the effective plastic strain rate. 

To compute the current value of the yield stresses needed for the yield function, in 
the developed model tabulated stress-strain curves are used to track the yield stress 
evolution. The user is required to input twelve stress versus plastic strain curves. 
Specifically, the required curves include uniaxial tension curves in each of the normal 
directions (1,2,3), uniaxial compression curves in each of the normal directions 
(1,2,3), shear stress-strain curves in each of the shear directions (1-2, 2-3 and 3-1), and 
45 degree off-axis tension curves in each of the 1-2, 2-3 and 3-1 planes. The 45 
degree curves are required in order to properly capture the stress interaction effects. 
By utilizing tabulated stress-strain curves to track the evolution of the deformation 
response, the experimental stress-strain response of the material can be captured 
exactly without any curve fit approximations. The required stress-strain data can be 
obtained either from actual experimental test results or by appropriate numerical 
experiments utilizing stand-alone codes. Currently, only static test data is considered. 
Future efforts will involve adding strain rate and temperature dependent effects to the 
computations. To track the evolution of the deformation response along each of the 
stress-strain curves, the effective plastic strain is chosen to be the tracking parameter. 
Using a numerical procedure based on the radial return method [7] in combination 
with an iterative approach, the effective plastic strain is computed for each time/load 
step. The stresses for each of the tabulated input curves corresponding to the current 
value of the effective plastic strain are then used to compute the yield function 
coefficients. 


DAMAGE MODEL OVERVIEW 


The deformation portion of the material model provides the majority of the 
capability of the model to simulate the nonlinear stress-strain response of the 
composite. However, in order to capture the nonlinear unloading and local softening 
of the stress-strain response often observed in composites [8], a complementary 
damage law is required. In the damage law formulation, strain equivalence is 
assumed, in which for every time step the total, elastic and plastic strains in the actual 
and effective stress spaces are the same [9]. The utilization of strain equivalence 


permits the plasticity and damage calculations to be uncoupled, as all of the plasticity 
computations can take place in the effective stress space. 
The first step in the development of the damage model is to relate the actual 
stresses to a set of effective stresses by use of a damage tensor M 
o=Mo,, (4) 
The effective stress rate tensor can be related to the total and plastic strain rate tensors 
by use of the standard elasto-plastic constitutive equation 
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where C is the standard elastic stiffness matrix and the actual total and plastic strain 
rate tensors are used due to the strain equivalence assumption. 

Given the usual assumption that the actual stress tensor and the effective stress 
tensor are symmetric, the actual stresses can be related to the effective stresses in the 
following manner, where the damage tensor M is assumed to have a maximum of 36 
independent components. 
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In order to maintain a one to one relationship between the effective stresses and the 
actual stresses (i.e. to ensure that a uniaxial load in the actual stress space does not 
result in a multi-axial load in the effective stress space), the damage tensor is assumed 
to be diagonal, leading to the following form. 
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However, an implication of a diagonal damage tensor is that loading the composite 
in a particular coordinate direction only leads to a stiffness reduction in the direction 
of the load due to the formation of matrix cracks perpendicular to the direction of the 
load. However, as discussed in detail in Goldberg, et al [5], in actual composites, 


particularly those with complex fiber architectures, a load in one coordinate direction 
can lead to stiffness reductions in multiple coordinate directions. To maintain a 
diagonal damage tensor while still allowing for the damage interaction in at least a 
semi-coupled sense, each term in the diagonal damage matrix should be a function of 
the plastic strains in each of the normal and shear coordinate directions, as follows for 
the example of the Mi; term for the plane stress case 
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Note that plastic strains are chosen as the “tracking parameter” due to the fact that, 
within the context of the developed formulation, the material nonlinearity during 
loading is simulated by use of a plasticity based model. The plastic strains therefore 
track the current state of load and deformation in the material. To explain this concept 
of damage coupling further, assume a load is applied in the 1 direction to an 
undamaged specimen. The undamaged modulus in the | direction is F,, and the 
undamaged modulus in the 2 direction isE,,. The stress-strain response of the 
material is assumed to become nonlinear (represented in the current model by the 
accumulation of plastic strain) and damage is assumed to occur. The original 
specimen is unloaded and reloaded elastically in the | direction. Due to the damage, 
the reloaded specimen has a reduced modulus in the 1 direction of E(''. The reduced 
area and modulus are a function of the damage induced by the loading and resulting 
nonlinear deformation in the | direction (reflected as plastic strain) as follows: 


Bae = (1 rae der )\E,, (9) 


where di; is the damage in the 1 direction due to a load in the 1 direction. 


Alternatively, if the damaged specimen was reloaded elastically in the 2 direction, due 
to the assumed damage coupling resulting from the load in the 1 direction, the 


reloaded specimen would have a reduced modulus in the 2 direction of E%j'. The 


reduced modulus is also a function of the damage induced by the load and resulting 
nonlinear deformation in the | direction as follows: 
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where d,, is the damage in the 2 direction due to a load in the 1 direction. Similar 
arguments can be made and equations developed for the situation where the original 
specimen is loaded in the 2 direction. 

For the case of multiaxial loading, the semi-coupled formulation needs to account 
for the fact that as the load is applied in a particular coordinate direction, the loads are 
acting on damaged areas due to the loads in the other coordinate directions, and the 
load in a particular direction is just adding to the damaged area. For example, if one 


loaded the material in the 2 direction first, the reduced modulus in the 1 direction 
would be equal to E“”’. If one would then subsequently load the material in the 1 
direction, the baseline modulus in the 1 direction would not be equal to the original 


modulus E,,, but instead the reduced modulus E¥*. Therefore, the loading in the 1 


direction would result in the following further reduction in the modulus in the 1 
direction: 
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These results suggest that the relation between the actual stress and the effective stress 
should be based on a multiplicative combination of the damage terms as opposed to an 


additive combination of the damage terms. For example, for the case of plane stress, 
the relation between the actual and effective stresses could be expressed as follows: 
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where for each of the damage terms the subscript indicates the direction of the load 
which initiates the particular increment of damage and the superscript indicates the 
direction in which the damage takes place. 


(11) 


CHARACTERIZATION OF DAMAGE MODEL 


There are two primary items needed for model characterization and input for the 
damage portion of the material model. First, the values of the various damage 


parameter terms d : need to be defined in a tabulated manner as a function of the 


effective plastic strain. Similar to the deformation model, the values of the damage 
parameters are defined in a tabulated, discretized form in order to reflect the actual 
material behavior in the most accurate manner possible. The values are tabulated as a 
function of the effective plastic strain in order to provide a unified framework to 
simultaneously track the evolution of multiple damage parameters under multiaxial 
loading conditions. Note that for the case of uniaxial loading the effective plastic 
strain equals the uniaxial plastic strain, which maintains consistency with the 
formulation described above. In addition to characterizing the damage parameters, the 
various input stress-strain curves need to be converted into plots of effective 
(undamaged) stress versus effective plastic strain in order to carry out the calculations 
required by the deformation (plasticity) model. As an example of how this process 
could be carried out, assume that a material is loaded unidirectionally in the 1 
direction. At multiple points, once the actual stress-strain curve has become nonlinear, 


the total strain (€,, ), actual stress (611), and average local, damaged modulus E41’ 
can be measured. Assuming that the original, undamaged modulus EF,, is known, 


since the damage in the 1 direction is assumed to be only due to load in the | direction 
(due to the uniaxial load), the damage parameters and effective stress in the | direction 
can be computed at a particular point along the stress-strain curve as follows: 


l-d}i= 11 
Ey, 
M,,=1-d,, 
ll Ul 
of — P11 (13) 
i > 
M,, 
eff 
p 07; 
an ll E 


These values need to be determined at multiple points, representing different values of 
plastic strain, in order to fully characterize the evolution of damage as the plastic strain 
increases. 

With this information, an effective stress versus plastic strain (<7) plot can be 
generated. From this plot, the effective plastic strain corresponding to the plastic 
strain in the | direction at any particular point can be determined by using the 
equations shown below, which are based on applying the principal of the equivalence 
of plastic work [7] in combination with Equation (2), simplifying the expressions for 
the case of unidirectional loading in the 1 direction [4]: 
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where ¢? is the effective plastic strain and de; is the increment of plastic strain in the 
1 direction. From this data, plots of the effective stress in the 1 direction versus the 
effective plastic strain as well as plots of the damage parameter d// as a function of 


the effective plastic strain can be generated. By measuring the damaged modulus in 


the other coordinate directions at each of the measured values of plastic strain in the 1 


direction, the value of the damage parameters d;,,d/;,d;, , etc. can be determined as 


a function of the plastic strain in the | direction, and thus as a function of the effective 
plastic strain. Ongoing efforts will involve developing and carrying out an appropriate 
experimental test matrix to characterize and validate the model for a series of 
representative composite materials. 


FAILURE MODEL 


A wide variety of failure models have been developed for composites. In models 
such as the Tsai-Wu failure model [6], a quadratic function of the macroscopic 
stresses is defined in which the coefficients of the failure function are related to the 
tensile, compressive and shear failure stresses in the various coordinate directions. 
This model, while mathematically simple and easy to implement numerically, assumes 
that the composite failure surface has an ellipsoidal (in 2D) or ovoid (in 3D) shape. In 
reality, composite failure surfaces often are not in the form of simple shapes. More 
complex models, such as the Hashin model [10], also utilize quadratic combinations of 


the macroscopic failure stresses, but utilize only selective terms in the quadratic 
function in order to link the macroscopic stresses to local failure modes such as fiber 
or matrix failure. However, an overall quadratic form to the failure functions (albeit in 
a piecewise fashion) are still assumed. This approach was extended in models such as 
those developed by Puck et al [11], Pinho et al [12] and Maimi et al [13], in which 
complex equations were developed to predict local failure mechanisms in terms of 
macroscopic level stresses. In this manner, the failure response and complex failure 
surfaces present in actual composites could be more accurately represented. However, 
in these advanced models, very complex tests are often required to characterize the 
model parameters and the applicability of the models may be limited to specific 
composite architectures with specific failure mechanisms. 

Within the context of transient dynamic analyses, failure models based on 
tabulated input are becoming more widely used. For example, in a model recently 
implemented in LS-DYNA for the impact analysis of metals [14], the effective plastic 
failure strain of the material is defined to be a function of the triaxiality, t, and the 
Lode parameter, 1, which are defined below: 
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In these expressions, J2 is the second invariant of the deviatoric stress tensor and J3 is 
the third invariant of the deviatoric stress tensor. The advantage of using the triaxiality 
and Lode parameter as the independent variables for defining the failure surface is that 
these terms can completely describe the load state of the material. Under plane stress 
loading conditions, the triaxiality can vary from a value of % for the case of biaxial 
tension to a value of -% for the case of uniaxial compression. Any other load 
conditions yield triaxiality values within these bounds. For example, the case of 
uniaxial tension yields a triaxiality value of 4, the case of uniaxial compression yields 
a triaxiality value of -'4, the case of pure shear yields a triaxiality value of 0 and the 
case of plane strain yields a triaxiality value of 1/V3 under plane stress conditions. 
The Lode parameter generalizes the process for the case of three dimensional loading 
conditions. The Lode parameter can vary between values of -1 and 1 depending on 
the state of stress that is applied to the material. Therefore, the load state of the 
material can be completely described by using these two parameters. By plotting the 
effective plastic strain of the material as a function of the triaxiality and the Lode 
parameter, the failure surface can be completely defined. These relationships are 
defined by using input curves and tables based on rigorously obtained experimental 
data as opposed to analytical functions, thus providing a tabulated representation of 
the failure surface. Advantages of this approach are that the experimental failure 
surface can be precisely defined in the model input and the variation in the material 
failure response due to varying load conditions can be accounted for in a systematic 
manner. Furthermore, the actual, experimentally obtained failure surface is used in the 
model, and a specific shape to the failure surface is not assumed. The failure model 


and relevant input are based on stress and strain invariants, which maintains the 
generality of the model. 

For the composite failure model developed in this work, the approach applied in 
the tabulated failure model for metals described above is adapted. The use of stress 
and strain invariants in composite failure models has been investigated by researchers 
such as Mayes and Hansen [15] and Feng [16]. Due to the fact that the failure 
response of composites can be brittle in certain material directions, instead of using the 
effective plastic strain as the dependent variable, a stress invariant first identified by 
Fleischer [17] is used as the dependent variable, defined as follows for the plane stress 
case (but generalizable for the full three dimensional case): 
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This invariant can be considered to be like a “radius” from the origin to the failure 
surface. The factor of 2 in front of the shear stress term reflects the symmetry of the 
stress tensor. This invariant can be plotted in terms of the triaxiality to define the 
failure surface under plane stress conditions. However, to account for the orthotropy 
of composite materials, the “radius” versus triaxiality curve needs to be defined for a 
variety of fiber orientation angles in order to fully describe the failure behavior of the 
composite. The angle can also be thought of as the angle between the first principal 
direction of the stress tensor and the structural axis system. An important point to note 
is that the failure surface (and the corresponding invariants) can be constructed on the 
basis of failure stresses or failure strains. To generalize to a three dimensional loading 
condition, the stress (or strain) invariant needs to be plotted as a function of the 
triaxiality and Lode parameter for a variety of fiber orientation angles. 

To characterize the failure model, the stress invariant defined in Equation (16) 
needs to be determined experimentally at the point of failure for a variety of loading 
conditions, expressed as appropriate values of the triaxiality and Lode parameter, for a 
variety of fiber orientation angles. These values need to then be expressed in a 
tabulated form. To utilize the failure model within a computational algorithm, given a 
computed stress state, the angle of the first principal stress direction and the baseline 
axis system, computed using the first eigenvector of the stress tensor, needs to be 
determined. Next, the triaxiality and Lode parameter of the stress condition need to be 
determined. Having then fully defined the “location” of the current load condition 
within the stress space, the stress invariant defined in Equation (16) needs to be 
computed. If the value of the invariant for the computed stresses for the selected 
angle, triaxiality and Lode parameter is less than the value at failure, failure is deemed 
to not have occurred. If the value of the invariant is equal to or greater than the 
corresponding failure value of the invariant for the particular load state, failure is 
considered to have occurred. 


CONCLUSIONS 


A generalized composite model suitable for use in polymer composite impact 
simulations has been developed. The model utilizes a plasticity based deformation 
model obtained by generalizing the Tsai-Wu failure criteria. A strain equivalent 
damage model has also been developed in which loading the material in a particular 


loading direction can lead to damage in multiple coordinate directions. A systematic 
approach to characterizing the parameters in the damage model has been developed. 
A methodology for utilizing a general, tabulated failure model has been formulated. 

Future efforts will included developing the detailed numerical algorithm to 
implement the damage and failure models into the material model being developed for 
inclusion within the LS-DYNA computer code. The failure model will also be refined 
to include appropriate techniques for element removal once the failure criteria has 
been satisfied. Methods to account for delamination failure will also be formulated. 
Extensive sets of verification and validation studies will also be undertaken in order to 
fully exercise the developed model. 
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